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A Lagrangian dynamic subgrid- 
scale model of turbulence 

By CL Meneveau 1 , T. S, Lund 2 AND W* Cabot 2 


A new formulation of the dynamic subgrid-scale model is tested in which the 
error associated with the Germano identity is minimized over flow pathlines rather 
than over directions of statistical homogeneity. This procedure allows the applica- 
tion of the dynamic model with averaging to flows in complex geometries that do 
not possess homogeneous directions. The characteristic Lagrangian time scale over 
which the averaging is performed is chosen such that the model is purely dissipa- 
tive, guaranteeing numerical stability when coupled with the Smagorinsky model. 
The formulation is tested successfully in forced and decaying isotropic turbulence 
and in fully developed and transitional channel flow. In homogeneous flows, the 
results are similar to those of the volume- averaged dynamic model, while in channel 
flow, the predictions are superior to those of the plane-averaged dynamic model. 
The relationship between the averaged terms in the model and vortical structures 
(worms) that appear in the LES is investigated. Computational overhead is kept 
small (about 10 % above the CPU requirements of the volume or plane- averaged dy- 
namic model) by using an approximate scheme to advance the Lagrangian tracking 
through first-order Euler time integration and lineax interpolation in space. 


1. Introduction 

The dynamic model (Germano et al 1991) for the parametrization of subgrid 
stresses in a Large-Eddy-Simulation (LES) is a means of utilizing information from 
the resolved turbulent velocity field fij(x, t) to dynamically compute model coeffi- 
cients. It is based on the algebraic identity, 

L{j = Tij — Tjj, (1) 

where 

Lij = UiUj — Tij = UiUj - and r t j = u75} — (2) 

Above, a tilde represents low-pass filtering with a filter- width of size A (comparable 
to the grid-size of the LES), while an overbar represents filtering at a scale 2 A. When 
the identity is written with the stresses Tij and Tij replaced by the Smagorinsky 
model, and Eq. (1) is enforced in a least-square error sense over all five independent 
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tensor elements (Lilly, 92), one obtains the following expression for the (dynamic) 
Smagorinsky coefficient: 


c»(x,0 = 


LpqMpq 

MpqMpg ’ 


( 3 ) 


where 

Mij = — 2A 2 [4|S|Sij - \S\S l} ] , (4) 


and Sij is the resolved rate-of-strain tensor. This version of the dynamic model in 
which the coefficient can vary from point to point is often referred to as the ‘local 
dynamic model’. 

There are two problems associated with the local dynamic model (Eq. (3)). First, 
as pointed out by Ghosal et al (1994), it is mathematically inconsistent to remove 
the coefficient from the filter operation (in r^J) as if it were a constant. Second, 
as observed during LES, during a-priori analysis of DNS data (Lund et a/., 1993) 
and when analyzing experimental data at high Reynolds numbers (Liu et a/., 1994; 
O’Neil & Meneveau, 1994), the coefficient field predicted by the local model varies 
strongly in space and contains a significant fraction of negative values. Negative 
values of c 2 are of particular concern because they lead to negative values of eddy 
viscosity in the Smagorinsky parameterization. This is destabilizing in a numerical 
simulation, and non-physical growth in the resolved velocity fluctuations is often 
observed (Lund et a/., 1993). 

Historically, the first problem was given very little attention while the second 
problem was dealt with by averaging terms in the equations for c 2 over space and/or 
time. When averaged, the numerator in Eq. (3) was generally found to be posi- 
tive, thus recovering the statistical notion of energy transfer to the subgrid scales. 
Averaging over homogeneous directions has been a popular choice, and excellent 
results were obtained in a variety of flows. As examples, Germano et al. (1991) 
and Piomelli (1993) average the equations over planes parallel to the walls in chan- 
nel flow simulations whereas Akselvoll and Moin (1993) average over the spanwise 
direction in a backward-facing step flow. While these averaging schemes proved to 
be effective at controlling possible instabilities and led to accurate results, rigorous 
justification for them was lacking. Additionally, homogeneity in either space or time 
was required. 

These problems as well as the lingering issue of extracting c\ from the filtering 
operation were addressed by Ghosal et al (1994) where a variational approach 
was used to account properly for the spatial variation of the coefficient within the 
filter operation. Using this approach, various prior models employing averaging 
were rigorously derived by imposing appropriate constraints in the solution to the 
variational problem. Finally, two stable local models were derived that did not 
make use of homogeneous directions. The first simply imposes the constraint that 
c l be non-negative. The second allows negative c 2 but enforces a budget for the 
reversed energy transfer through inclusion of a subgrid-scale kinetic energy equation. 
These latter two models have been tested in a variety of flows and are applicable to 
complex geometry flows under unsteady conditions. 
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While the work of Ghosal et al. (1994) has provided rigorously- derived methods 
applicable to inhomogeneous flows, there is still room for improvement. The con- 
straint c 2 s > 0 is hard to justify on other than heuristic grounds, and the numerical 
solution of the integral equation can be expensive. The kinetic energy formulation 
removes the conceptual problem associated with the constraint c 2 > 0, but only at 
the additional expense of two more integral equations and one transport equation. 
Also, new constraints for model coefficients in the kinetic energy equation have to be 
introduced. Therefore, schemes that make use of averaging continue to have appeal 
due to their demonstrated accuracy and relative ease of implementation. At the 
same time, current averaging schemes require at legist one homogeneous direction 
(in space or time) and are therefore not applicable in fully inhomogeneous unsteady 
flows. If this restriction could be removed, an equally general and perhaps simpler 
alternative to the integral equation of Ghosal et al. (1994) would be available for 
inhomogeneous flow simulations. 


The objective of this work is to develop a simple, but generally applicable, aver- 
aging scheme. As originally suggested by O’Neil & Meneveau (1993), we propose 
to average over particle trajectories rather than directions of statistical homogene- 
ity. Particle trajectories are always well defined objects that in no way rely on 
special boundary conditions or assumptions of statistical homogeneity. Since parti- 
cle trajectories are the natural directions associated with fluid flow, averaging the 
equations for c 2 s over these directions has some physical appeal. It is reasonable to 
expect that turbulent eddies evolve along particle paths and that the turbulence 
energy cascade should be most apparent when viewed in a Lagrangian coordinate 
system. As reported by Meneveau & Lund (1994), there is evidence to suggest that 
this is indeed the case. If the energy cascade does in fact proceed mainly along fluid 
trajectories, then it would seem logical to postulate that the subgrid- scale model 
coefficient at a given point x should depend in some way on the history of the flow 
along the trajectory leading to x. This picture should be contrasted with that of 
conventional schemes where spatial averaging removes the local details of the flow 
structure and the turbulence development history is completely ignored. Eulerian 
time averaging suffers from similar deficiencies since the advection of structures is 
ignored. 


The Lagrangian model is derived by requiring that the error in Germano’s identity 
be minimized along fluid trajectories. This procedure leads to a pair of relaxation 
transport equations that carry the statistics forward in Lagrangian time. We show 
that these equations can be solved in an approximate fashion in a numerically 
efficient way. The model is applied to a variety of test cases including forced and 
decaying isotropic turbulence, fully developed channel flow, and transitional channel 
flow. In each case, the model is shown to produce results equal or superior to those of 
spatially-averaged versions of the dynamic model. At the same time, the numerical 
solutions to the transport equations increase the computational workload by only 
about 10% as compared with the spatially- averaged approach. 
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2. The Lagrangian dynamic model 

2.1 Formulation 

We propose to determine the model coefficient c^(x, £) by minimizing the error 
in Germano’s identity along particle trajectories. Consider a particle located at 
position x at time t. The trajectory of this particle for earlier times t* < t is 


i(t) = x — u [z(t ff ),t n ]dt n 


(5) 


The error associated with Germano’s identity at any point along the trajectory is 

e«j(M') = (6) 

Here we have assumed that c 2 varies negligibly in space over the scale of the test 
filter and have therefore removed it from the filter operation. 

The total error is defined as the pathline accumulation of the local error squared, 


i 

E = J e 0 (z (<'),0 W(t - t') dt' 


(7) 


The weighting function W(t — t 9 ) is introduced here in order to control the relative 
importance of events near time t with those of earlier times. As described below, 
we shall weight the error at time t most strongly and assign a decreasing weight to 
earlier times. It is now assumed that c 2 s varies negligibly in time over the scale of 
the weighting function W(t — t'). In this case, c*(z,t') may be replaced by c^(x,£), 
and the total error is then minimized with respect to c 2 by enforcing 


dE r 

= J 


dt 




del 


W(t - t') dt' = 0. 


( 8 ) 


Making use of Eq. (6) (with cl(z,t') replaced by c^(x, <)) and solving for cj, one 
then obtains 

2, i \ 7 lm 

c;(M) = ^ — , 

-l-MM 


where 


t 

Tlm(x,<) = J LijMij(z(t'), t') W(t - 1') dt', 


— OO 

t 


(9) 


( 10 ) 


Imm(x,<) = J MijMij(z(t'),t') W{t — t') dt', (11) 
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FIGURE 1 . Sketch of fluid trajectory of the resolved LES velocity field. The error 
associated with the Germano identity is weighted with an exponentially decreasing 
function (indicated as different gray levels), backwards in time, to yield a current 
value for the model coefficient at point x and time t. 


The function W(t — t') is a free parameter in the current formulation, essentially 
defining the extent backward along the pathline over which we choose to minimize 
the error. Although several appropriate weighting functions are possible, an expo- 
nential weighting of the form W(t - t') — j— 1 e -( t -‘')/u has the distinct practical 
advantage that the integrals Tim and Tmm are solutions to the following relaxation 
equations 


DTj m 

Dt 

DTmm 

Dt 


dlLM , - X7T 

~dT + uVIlm 

MM 


— {LijMij -Ilm), 
^ ( MijMij — Tmm) ■ 


( 12 ) 

(13) 


In the context of LES, solving such transport equations is much more natural than 
having to perform integrals backwards in time according to Eqs. (10) and (11). 
Fig. 1 illustrates the basic idea of averaging over pathlines with an exponentially- 
decreasing memory. 


2.2 Relaxation time scale 

The time-scale T controls the memory length of the Lagrangian averaging, and 
several choices can be made. The model coefficient should be responsive to changes 
that occur on the time-scales associated with the smallest resolved turbulent mo- 
tions. Thus, one could choose T based on variables at the grid-scale. Some possible 
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choices are the following: (a) T ~ (b) T ~ |S|-\ (c) T ~ A 

(d) T ~ A(L ii M’ 0 -)" 1/4 . ( e ) r ~ AJ" 1 //, and (f) T ~ AT"]/ 4 . In fully devel- 
oped turbulence, all of these time-scales are of the same order of magnitude on 
average. The first four choices are based on local values, which would mean that 
T is a strongly fluctuating variable while (e) and (f) are based on the smoother, 
Lagrangian averages themselves. 

Option (f) has several attractive features. Physically, it can be interpreted as 
a time scale for energy flux since it is formed by contracting a stress L with the 
strain-rate like tensor M. Thus, it may be indicative of the speed at which energy is 
being cascaded towards the grid-scale. Furthermore, if LijMij < 0 for a persistent 
time along the pathline, then Tim approaches zero. T evaluated according to (f) 
then tends to oo, i.e. the memory time increases. In other words, the current 
values are weighted less and less strongly relative to the past ones, if they are of the 
backscattering type. This is useful in the implementation with the Smagorinsky 
model where we wish to restrict the Smagorinsky expression for the modeling of 
energy dissipation only. The Germano identity is thus weighted much less heavily 
when LijM t j < 0 in a persistent fashion, i.e. we opt for ‘learning’ as little as possible 
about the coefficient from the resolved field when it would predict backscatter. 

Equation (12) can now be written as 

= |lL (14) 

where 9 is a dimensionless coefficient of order unity. If Tim reaches zero, its rate 
of change is zero as well. Therefore, Tim cannot become negative, and the result- 
ing dynamic model will not suffer from numerical instability due to negative-eddy 
viscosities. We point out, however, that if LijMij < 0, the approach of Tlm to 
zero is not exponential, but of the power-law type (as (to — t) 4 / 3 ). This means that 
after the (finite) time at which Tim = 0, the solution becomes complex. Thus, 
in practice, the solution must still be ‘clipped’ to zero during such times. This type 
of clipping is much less drastic than previous approaches since the coefficient field 
approaches zero smoothly with zero slope. 

A judicious choice for the dimensionless coefficient 9 must now be made. In- 
tuitively, we must average over a few ‘events’ of the variable LijMij along the 
pathline. The average duration of such events is expected to be of the order of 
A < LijMij >“i, but in order to quantify this assertion, we analyze results from 
DNS of forced isotropic turbulence. The goal is to compute the Lagrangian auto- 
correlation function of the scalar variable LijMij. The method employed to follow 
fluid trajectories and to compute the autocorrelation is the same as described in 
Meneveau & Lund (1994). For comparison, we also compute the Lagrangian auto- 
correlation function of the scalar |M| 2 as well as their Eulerian fixed-point two-time 
autocorrelation functions. The Lagrangian tracking was done in a sequence of DNS 
velocity fields computed on a 128 3 mesh. The ensemble has a microscale Reynolds 
number of R\ = 95.8. Each field was filtered with a Fourier cutoff filter at a scale 
corresponding to 4 mesh spacings. Lagrangian and Eulerian autocorrelations were 
then computed for quantities derived from the filtered velocity fields. 
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FIGURE 2. Lagrangian and Eulerian autocorrelation functions calculated from 
a filtered DNS of forced isotropic turbulence. • : Lagrangian autocorrelations of 

LijMij ; *• : Lagrangian autocorrelations of MijMij ; : Eulerian temporal 

autocorrelations of LijMij ; : Eulerian temporal autocorrelations of MijMij . 

Fig. 2 shows the computed autocorrelations. Time is non-dimensionalized based 
on the space-averaged value of LijMij . As expected, the Lagrangian autocorre- 
lations decay at a slower rate than the Eulerian ones, but the difference is small 
due to the fact that the mean velocity of this flow is zero. Also, the decay of the 
LM and MM terms is quite similar. The main observation is that after a time 
~ 2A < LijMij >”S the autocorrelation almost vanishes. This suggests that av- 
eraging over Lagrangian time spans equal to this interval is sufficient to smooth 
instantaneous fluctuations. In summary, during the present work we choose 

T = 2AX‘* (15) 

as the time-scale characterizing the exponential memory with which the Germano 
identity is enforced. 

2.3 Numerical method 

In principle, the implementation of the Lagrangian dynamic model requires the 
solution of two additional transport equations (Eqs. (12) and (13)) during the LES. 
This undoubtedly increases the computational cost associated with the subgrid 
modeling. However, the considerable flexibility of choice of the averaging domain 
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suggests that high numerical accuracy in solving Eqs. (12) and (13) may be unnec- 
essary. Therefore, we use a particularly simple formulation based on discretizing 
Eq. (12) in time as follows: 



Eq. (13) is dealt with in a similar manner. Positions x are coincident with grid 
points of the simulation. The value of X"^ at the previous time-step and at the 
upstream location x — u n A t can be obtained by multilinear interpolation. Finally, 
the new values at the grid points are solved for. The result is a weighted sum of 
the interpolated prior value and the current source term at the grid point: 

I££(x) = H { 6 [LyJIfyr+^x) + (1 - 6)I£ m (x - u"At) } (17) 


and 

I?m( x ) = '[MijMij ] n+1 (x) + (1 - e)IJf M (x - u n A<), 

where 

e ~ TTtefF ' and r n = 2A(i£ M r*, (is) 

where H{x} is the ramp function ( H{x } = x if x > 0, and zero otherwise). The 
ramp function is introduced to clip the solution away from complex values. 

Finally, we point out that the process of spatial interpolation between grid points 
introduces some numerical diffusion to the fields Ilm and Imm* Physically, such 
diffusion effectively 'thickens’ the pathline over which the averaging is being per- 
formed, but this would not seem to be a worrisome aspect for this model. 

2.4 Statistical features of the model 

As a next step, the model is implemented in a LES for the simulation of forced 
isotropic turbulence on a 32 3 grid. The code is a variant of the pseudo- spect ral 
method developed by Rogallo (1981). Forcing is achieved by holding the Fourier 
amplitudes fixed within the sphere A: <2. Test filtering is achieved through a 
Fourier cutoff at twice the grid scale. 

The velocity field is initialized in the usual manner by superposing Fourier modes 
with a prescribed spectrum but random phases, and projection onto the divergence- 
free space. Additionally, initial condition for the fields Jlm and Imm must be 
prescribed. For initializations corresponding to turbulent flows, we propose to set 


Jmm(x, 0) = MijMij(x,0), Ilm(x, 0) = c^(0) M o M o (x,0), (19) 

where c*(0) = 0.16 is the traditional value of the Smagorinsky constant. Thus at 
the initial time, the model involves a position-independent, prescribed coefficient. 
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FIGURE 3. Probability density functions of the coefficient c 2 , computed from the 
Lagrangian dynamic model in a pseudo- spectral LES of forced isotropic turbulence 
at Re = oo, on a 32 3 mesh. The circles represent the asymptotic pdf long after 
initial transients have died out (obtained here at t = 6 < T >, where < T >= 
2A < LijMij Evolving pdfs on both sides illustrate delta function pdfs 

with the wrong initial conditions quickly reach the asymptotic statistics. Curves 
that peak to the left of the asymptotic curve correspond to c*(0) = 0.005; those 

peaking to the right evolve starting from c\ = 0.075. : (1 time-step); : 

(t = 0.1 < T >); : (t = 0.45 < T >); : (t = 0.95 < T >); 

(t = 1.9 < T >). For reference, in this simulation the time scale associated with 
the resolved strain-rates was < SijSij >~ l / 2 = 0.26 < T >. 


For initializations corresponding to laminar flows, we propose to set c 9 — 0 in the 
above expressions. 

When the LES of forced isotropic turbulence is started, fluctuations of the La- 
grangian dynamic coefficient c 9 quickly build as different values of LijMij begin to 
affect the averages. Once a statistical steady-state has been reached, these fluctua- 
tions are characterized by the probability density function of the coefficient shown 
by solid circles in Fig. 3. Notice the small spike at c 9 = 0, arising from the regions 
in which c a is clipped at zero, away from complex values (on about 5% of the points 
in this case). Initial transients leading to such a steady-state distribution are rela- 
tively short. This can be appreciated by observing the time development of the pdfs 
when the ‘wrong’ initial condition is employed for c*(0). In one case, we start with 
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LijM t j, Iim 

FIGURE 4. Probability density functions of numerators. : evaluated locally 

( LtjMij ); : after Lagrangian averaging ( Ilm ). These distributions are cal- 

culated from a 32 3 -node, pseudo- spectral LES of forced isotropic turbulence that 
uses the Lagrangian dynamic Smagorinsky model. To increase the sample, pdfs are 
accumulated over several independent fields. 




MijMij , 1mm 


FIGURE 5. Probability density functions of denominators. : evaluated 

locally ( ); : after Lagrangian averaging {1m m)- Details as in Fig. 4. 
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FIGURE 6. Probability density functions of model coefficients taken from a 32 3 

Lagrangian dynamic model LES of forced isotropic turbulence. : coefficient 

evaluated locally; : coefficient from the Lagrangian model. 


c^(0) = 0.005, and in another case with c^(0) = 0.075. In both cases the asymptotic 
distribution is reached after times of the order of 2 < T > where T is the time scale 
defined by Eq. (15). We conclude that the proposed method of initialization is 
acceptable since the simulation ‘forgets’ the initial state after only few grid-scale 
turnover times. This is comparable to the time it takes the simulation to build up 
realistic phases in the resolved velocity field, starting from the random-phase initial 
condition. To further document the effect of the Lagrangian averaging, we com- 
pute the probability density functions of Tim and Tmm and compare them with 
those of the local values LijMij and Figs. 4 and 5 show these results. As 

expected, the distributions become narrower after the Lagrangian averaging. By 
construction, there are no negative values of Tim (Fig. 4). In terms of denomina- 
tors, the averaging is seen to virtually eliminate values near zero. The pdf of Tmm 
approaches the origin with negligible slope while the probability of the local value 
of MijMij being close to zero is considerable. In Fig. 6, we show the measured 
pdf of the coefficient c 2 s itself. As can be seen, the variance of the coefficient in the 
local formulation is greatly reduced by the Lagrangian averaging. Also, no negative 
values exist although a finite number of points 5%) exhibit c 2 = 0 as indicated 
by the spike at the origin. 
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2.5 Relationship to vortical flow- structures 

The goal of this section is to make qualitative observations pertaining to possible 
relations between the terms Ilm, 1mm , and discernible flow structures that may 
appear in the resolved velocity field during LES. 

First, we report the existence of tubular structures that characterize regions in 
which the resolved vorticity vector has high magnitude, in our 32 3 LES. Fig. 7a 
shows iso-surfaces of vorticity magnitude (at a threshold of |u;|^ = 2.4 < u; 2 > 1 / 2 ). 
Clearly, Tat worms 7 exist in the solution. The existence of tubular vortical structures 
in LES has also been observed recently by Briscolini and Santangelo (1994), using 
a different subgrid model. One interesting question to be answered is whether the 
prediction of such Tat worms 7 by LES is realistic. We recall that DNS predicts 
worms with very small diameters of about four Kolmogorov scales (Jimenez et al . , 
1993)). Surely they cannot be captured by a LES at Re = oo. The relevant question 
is whether a field generated by DNS and then low-pass filtered at inertial-range 
scales comparable to the LES grid-size exhibits Tat worms 7 that are comparable 
to those predicted by LES. We have performed such an operation based on the 
128 3 forced DNS described earlier and have visualized regions of high vorticity- 
magnitude. We indeed observed Tat worms’ that were of similar appearance than 
those of the LES (see also Fig. 17 of Vincent & Meneguzzi, 1991). It must be 
recognized that the ‘high-vorticity 7 regions in the filtered DNS correspond to much 
lower vorticity-magnitudes than those of the unfiltered fields. This is the reason 
why these Tat worms’ are not visible when analyzing the unfiltered DNS fields. 
In summary, we observe elongated vortical regions in LES and believe that their 
existence is a realistic prediction by the simulation since they also exist in low-pass 
filtered DNS fields. 

The next issue to be addressed is whether the Lagrange- averaged quantities that 
enter our dynamic model bear any relationship to such local structures. Fig. 7(a) 
shows contour plots of Ilm on two planes of the computational cube, selected to 
cut some of the most visible vortical structures. The field is chosen at some time 
long after the simulation has reached statistical steady-state (15,000 time-steps). 
Fig. 7(b) shows a similar graph for 1mm- Generally, it is apparent the contours of 
both Ilm and 1mm are somewhat ‘correlated’ with the presence of worms. The 
contours peak in the neighborhood of the worms while not much activity is seen 
in regions that are far removed from the structures. Upon closer examination, we 
observe that the peaks in Ilm and 1mm are most often located near the cores of 
the worms but not inside of them. Many times the maximum values occur between 
two closely spaced worms. These are expected to be regions of large straining 
and turbulence generation. Also, considerable correlation is seen between both 
fields Ilm and 1mm (It is p{1lm ,1mm) ~ 0.8 while the correlation between the 
local values is much lower, p(LijMij, MijMij) ~ 0.4). This increased correlation 
between numerator and denominator is instrumental in decreasing the variance of 
the predicted model coefficient. 

Clearly, a detailed understanding of the relationship between the coefficient c 2 
and local flow structures, and of their dynamical interplay and relevance, is still 




FIGURE 7. Visualization of high-magnitude vorticity regions in LES of isotropic 
forced turbulence at Re = oc. Surfaces correspond to points at which \u\ = 2.4 < 
cl? 2 > 1 ! 2 . On planes, contours of different variables shown, (a) Tim, (b) Imm? 
(c) \S\ (proportional to eddy viscosity with volume-averaged coefficient) and (d) 
(Jlm /Zmm) |S| (proportional to eddy viscosity computed from Lagrangian dy- 
namic model). 
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elusive. Nevertheless, we have shown that the Lagrangian averaging preserves some 
spatial locality in the model. Spatially localized events in the numerators and 
denominators used to compute the model coefficient bear some relationship to lo- 
cal flow structures. The volume- averaged dynamic model would have generated a 
position-independent coefficient that is oblivious to local flow structures. To observe 
the effect on the predicted eddy viscosity, Figs. 7(c) and 7(d) show contour plots 
of the strain-rate magnitude \S\ and of the expression (Xlm/^A/M )|* 5|* The former 
is proportional to the eddy viscosity predicted with a volume averaged coefficient 
while the latter is proportional to the eddy viscosity predicted by the Lagrangian 
dynamic model. Both show peaks surrounding the worms, but the precise location 
of these peaks differs. Also, the Lagrangian dynamic eddy viscosity appears to be 
more intermittent in the sense that more eddy viscosity is concentrated near the 
structures while being lower and fluctuating less in the regions far away from the 
structures. 

3. Applications 

In this section, we report applications of the Lagrangian dynamic model to several 
test-cases. We consider forced and decaying isotropic turbulence and channel flow. 
These flows could have also been treated with the traditional dynamic model with 
averaging over statistically homogeneous directions (alid they have in the past). Our 
purpose in choosing these simple flows is to test the model in well-understood cases 
and show that good results can be obtained. This is a necessary first step before 
applications to unsteady and complex- geometry flows should be attempted where 
many other effects such as numerics, etc. may influence the results and obscure the 
role of the subgrid model. 

3 A Forced isotropic turbulence 

LES of forced isotropic turbulence is performed on both 32 3 and 128 3 grids, using 
the code already described in section 2.4. The simulation is run for 15,000 and 6,600 
time-steps on the 32 3 and 128 3 grids, respectively. 

Figs. 8(a) and 8(b) show the resulting radial energy spectra. The wavenumbers 
and energy density are made dimensionless with the grid wavenumber and the aver- 
aged subgrid-scale energy dissipation (— < SijTij >). Figure 8(b) is premultiplied 
by & 5 / 3 . In these 4 mesh- Kolmogorov units’, one expects simulations with different 
meshes to collapse at high wavenumbers, and the spectra to follow the universal 
power- law in the inert ial range. The dotted line in Fig. 8(a) shows a power-law 
(k/k A y*/\ A slight decay below the power law for k/k& > 0.3 and a ‘pile-up’ very 
close to the cut-off wavenumber k& are visible. These are known effects of physical 
space eddy viscosity closures, which do not have a ‘cusp’ near These defects 
appear not to be remedied by the dynamic model in its Lagrangian implementa- 
tion. We have confirmed that the same is true for the traditional dynamic model 
by running the same program with the volume- averaged dynamic coefficient. 

With regard to computational cost, we find that the CPU time for the simulation 
with the Lagrangian averaging was higher by about 9% when compared to that of 
the volume- averaged dynamic model. Most of the additional time was spent in the 
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(°) 



(b) 



FIGURE 8. Radial energy spectra of LES using Lagrangian dynamic Smagorinsky 
model at Re = oo. Wavenumbers are non-dimensionalized with grid-wavenumber 
while non-dimensionalization of energy density also involves mean subgrid-scale 
energy dissipation rate, e = — < TijSij >, computed from simulation, (a) conven- 
tional spectra, (b) premultiplied spectra. : 128 3 simulation; : 32 3 

run; : slope —5/3. 32 3 spectrum averaged over 374 independent samples 

taken from 15,000 time- steps (approximately 13 integral large-scale turnover times 
Lju\ with L = 27r). The 128 3 spectrum is based on 290 samples taken from 6,600 
time-steps (approximately 4 large-scale turnover times). 





286 


C. Meneveau et al. 



Figure 9. Temporal decay of turbulent kinetic energy in isotropic turbulence. 

: 32 3 LES using the Lagrangian dynamic model; • : filtered experimental 

results (Comte-Bellot &; Corrsin, 1971) in grid turbulence which decays downstream 
of a grid. C/q and M are the mean fluid speed and the spacing of the turbulence- 
generating grid in the experiment. 

linear interpolations. Two additional scalar arrays had to be defined, for Tlm and 
2 mm* Compared to overall memory requirements, this addition was not significant. 

8.2 Decaying isotropic turbulence 

In order to test the model in an unsteady case, we perform LES of decaying 
isotropic turbulence. Meaningful comparisons can then be made with the experi- 
mental results of Comte-Bellot &: Corrsin (1971). The initial 3-D energy spectrum 
is made to match the experimental measurements at their earliest time. The phase 
of the Fourier coefficients is chosen to be random so that the initial velocity field 
had Gaussian statistics. The dimension of the computational box is chosen to be 
roughly 4 integral scales. 

Fig. 9 shows the decay of the kinetic energy compared with the experimental 
results of Comte-Bellot &; Corrsin (1971). The predicted initial decay appears to 
be a little slower than the experimental rate, but the overall agreement is good. 
Of course the agreement could have been improved by using a slightly larger value 
of c^(0) as initial condition - an after-the-fact adjustment that we opted to avoid. 
A comparison of the spectra at the three different times at which experimental 
results are available is shown in Fig. 10. The decrease in overall kinetic energy and 
the decrease of k at which the spectra peak (increasing integral scale) is clearly 
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FIGURE 10. Radial energy spectra for decaying isotropic turbulence at different 

times. : 32 3 LES with Lagrangian dynamic model; • : experimental results 

of Comte-Bellot & Corrsin (1971). Scaling parameters are defined in Fig. 9. L — 
10. 8M is the computational box size. 


reproduced well. 

We conclude that the model is able to reproduce important features of this time- 
dependent flow. 

3.3 Fully developed channel flow 

In this section we describe the application of the Lagrangian dynamic model to a 
pseudo-spectral simulation of plane channel flow. For comparison, another LES is 
performed with the traditional implementation of the dynamic model in which the 
terms axe averaged over planes parallel to the wall. The flow Reynolds number is 
selected to match the experimental data by Hussain & Reynolds (1970) to permit 
detailed comparison. 

The channel flow simulations are performed with in a pseudo-spectral code (Kim 
ei a/., 1987) in a numerical domain with streamwise, wall-normal, and spanwise 
dimensions of 37r x 2 x 37r/4 (in units of channel half- width d) on a 48 x 65 x 64 
mesh. Chebyshev polynomials are used in the wall- normal direction on a collocated 
grid; Fourier transforms are used in the homogeneous streamwise and spanwise 
directions on a uniform grid. Real space (tophat) filtering is used for the dynamic 
test filtering procedure and is performed explicitly only in horizontal planes. The 
equivalent filter width A eq is taken to be the geometric mean of unidirectional grid 
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FIGURE 11. Mean velocity profiles in fully developed channel flow. : 

LES with plane- averaged dynamic model; : LES with Lagrangian dynamic 

model; • : Experimental measurements of Hussain & Reynolds (1970). (a) Wall 
distance in units of the half- channel width d, (b) in wall units. 
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FIGURE 12. Profiles of second-order moments of the resolved velocity field. : 

from the plane- averaged LES; P : u t rm3 from Lagrangian model; o : u' rmj 
from the experimental measurements of Hussain & Reynolds (1970). o : Spanwise 

w* rm9 from Lagrangian LES; : from plane-averaged model. A : Wall-normal 

v' rrns from Lagrangian LES; : from plane- averaged model. ★ : Resolved 

shear stress < uV > from Lagrangian LES; : from plane-averaged model. 

widths (this procedure is justified for moderate grid anisotropies as shown in Scotti 
et <z/., 1993). For the plane- averaged LES, averaging of the dynamic coefficient is 
performed in horizontal planes. 

The approximate Lagrangian interpolation for the horizontal directions is im- 
plemented in this code as described in §2.3. The wall-normal direction requires 
different treatment due to the stretched mesh used in that direction. The trans- 
formation 8 = cos ~ l (y/d) is used to map the stretched mesh into a uniform one. 
The wall-normal advection term, vd/dy , is recast as v$d/d8 and the interpolation 
is performed in 8 identically to the horizontal directions, but using v $ = — vj sin#. 
The wall planes are treated specially with Tim = 0 and 1mm approximated by 
values at the nearest off- wall plane. There was also the possibility that the interpo- 
lation might attempt to place approximated points at the previous time step beyond 
the walls; however, the CFL condition gives sufficiently small time steps that this 
situation is never encountered. 

A target friction Reynolds number Re r (~ u r d/i^, where the friction speed u T 
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is the square root of the mean total wall stress, and v is the molecular viscosity) 
of 650 was chosen, corresponding to one set of experimental data by Hussain &; 
Reynolds (1970). The channel flow is started from a flow field at lower Reynolds 
number and is allowed to evolve to near statistical equilibrium, with Re r « 641 
in the last runs. The initial conditions for Jim and Xmm are chosen as in the 
homogeneous case but with c^(0) as function of y matching the values of a previous 
plane- averaged dynamic simulation. Using the Lagrangian formulation proved to 
be more expensive than the standard plane- averaged method by 10% in CPU due, 
in part, for the need to perform a division at each point to compute the dynamic 
coefficient rather than at each plane. The Lagrangian method also requires extra 
mass storage of Jim and Jmm between runs. 

The averaged statistics will be shown first, followed by a more detailed analysis 
of additional variables. 

Fig. 11 shows the mean velocity profile in the half-channel, in outer units (a) and 
wall units (b). As can be seen, at the resolution of the present LES, the plane- 
averaged model predicts an excessive center-line velocity (smaller losses) for the 
prescribed pressure gradient. The Lagrangian model yields a centerline velocity 
slightly below the measured values although the magnitude of the error is consider- 
ably smaller than that of the plane-averaged case (6.8% error in centerline velocity 
for the plane- averaged model and -1.8% for the Lagrangian model). Fig. 12 shows 
the profiles of rms velocities and Reynolds shear stress of the resolved fields, and 
a comparison of the rms streamwise velocity with the measurements of Hussain & 
Reynolds (1970). In the core region (for y/d > 0.2), the LES with both models 
fall below the measured values to a large extent because the former do not include 
the subgrid portion of the energy. Closer to the wall, both LES erroneously over- 
predict u f rms , but the Lagrangian model does a better job than the plane-averaged 
one. Interestingly, the magnitude of the resolved shear stress for the Lagrangian 
model is larger than that of the plane- averaged case. This is possibly the cause for 
the increased (more realistic) losses in the Lagrangian simulation. The mean eddy 
viscosity predicted by both LES is shown in Fig. 13. It is computed according to 

< Vt > (y) =< c](x,y,z,t ) 2 A e ,|5| > IiZ , (20) 

where the averaging is performed along x y z planes and over several times. The 
coefficient c\ is either computed according to the plane-averaged or the Lagrangian 
dynamic model. It can be seen that over much of the log-layer, the Lagrangian 
model generates a lower eddy viscosity compared to the plane- averaged dynamic 
model. We have checked that this reduction is due primarily to a decrease in the 
dynamic coefficient c 2 as opposed to reduced strain-rate magnitudes. The reduced 
eddy viscosity is likely to generate less SGS dissipation of resolved turbulence, which 
in turn is probably the cause for the increased resolved shear-stress observed before. 

As an aside, an important feature of the dynamic model is that it exhibits the 
proper near- wall scaling for the SGS eddy viscosity when the sublayer is numerically 
resolved (Germano et a/., 1991), namely u t ~ (y + ) 3 . As can be seen in Fig. 13, 
this scaling is followed quite well by the plane- averaged case (as observed before by 
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FIGURE 13. Ratio of mean eddy viscosity to molecular viscosity taken from 
the channel flow simulation, o : Lagrangian dynamic model; □ : plane- averaged 
dynamic model; : v% ~ (y + ) 3 power-law. 


Germano et al 1991; Piomelli, 1993). The mean eddy viscosity from the Lagrangian 
model also decays very quickly but at a somewhat slower rate (approximately as 
v t ~ (y + ) 2,5 in our case). We shall return to this issue at a later stage. But we 
stress that near the wall such minute differences are unlikely to have any practical 
effect since there the molecular viscosity strongly dominates. With the purpose 
of documenting the statistics of the model coefficient c 2 8 and its evolution away 
from the initial condition, we show in Figs. 14(a)-(c) probability-density-functions 
of c] at different times and different elevations from the wall. The pdf at t = 0 
is a delta - function at the plane- averaged value of the dynamic coefficient, which 
is used as an initial condition. As can be seen for t/ + = 641 and = 12 (core 
and near- wall region), the convergence of the pdf to the asymptotic value (circles) 
is nearly complete after 80-160 time- steps. This duration corresponds to about 
v/ul rss 25 — 50 viscous times or dju T ~ 0.04 — 0.08 outer times. At — 108, 
the convergence is slower because in the log-layer the initial guess for c 2 a is worst. 
Still, after between 300 and 600 time-steps, the asymptotic state is reached for 
the fluctuations in c 2 s . Figs. 14(b) and (c) clearly show the considerable decrease 
of typical c\ values in the Lagrangian model as compared to the plane- averaged 
model. 
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FIGURE 14. Probability-density functions of dynamic model coefficient cj for 

different time steps, n. : n=2; : n=10; : n=20; : n=80; 

: n=320; : n=640; o : n ~ 38,000. Different graphs are at different 

heights above the wall: (a) is at = 641, (b) at = 108 and (c) at y + = 12. 
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LES - plane- averaged 

y + 

< LijMij > 

< M^Mij > 

12 

2.912 10 4 
3.054 10 7 

108 

5.362 10 2 
4.932 10 4 

640 

1.099 10 1 
6.998 10 2 

LES - local Lagrangian 

y + 

< ?lm > 

< %MM > 

12 

2.690 10 4 
2.814 10 7 

108 

7.076 10 2 
2.537 10 5 

641 

1.056 10 1 
7.777 10 2 


TABLE 1. Numerators and denominators in the expressions for the dynamic 
coefficients averaged over sample planes. 



FIGURE 15. Scatter plot of 1mm versus vertical velocity v 9 in LES of channel 
flow, using the Lagrangian dynamic model. 

The main issue left to answer is why the Lagrangian model generates such de- 
creased coefficients in the log-layer. For this purpose, the average values of numer- 
ators and denominators are evaluated separately for both models on some sample 
planes as given in Table 1. 

The largest discrepancy can be seen by comparing the denominators < MijMij > 
and < 1mm > at = 108. A possible reason for this discrepancy can be deduced 
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by comparing < > at = 12 and at y + = 108. < MijMij > and 

< T'MM > a re several orders of magnitude higher in the near- wall region, as is to 
be expected for a variable based on the strain-rate (to the fourth power). During 
ejection events, fluid particles that were close to the wall reach deep into the log- 
layer, thus convecting elevated values of 2mm upwards. This feature can be deduced 
from Fig. 15, which shows a scatter plot of the Lagrangian denominator 2 mm as 
function of the local vertical velocity. Clearly, large values of 2m m are associated 
with positive values of v* , which are indicative of ejection events or bursts. 

The net effect is that the Lagrangian model is less dissipative as far as bursts 
are concerned. They can survive longer and feed more turbulence into the channel 
flow, producing more realistic (higher) levels of Reynolds-averaged (resolved) eddy 
viscosity and losses. 

It is likely that a similar phenomenon causes the near-wall scaling of Lagrangian 
eddy viscosity to be less steep than that of the plane- averaged model. Occasionally, 
‘sweeps’ bring log-layer material into the sublayer and effectively increase the model 
coefficient and eddy viscosity above that of the plane-averaged model. Numerical 
diffusion is also likely to play a role in reducing spatial differences in 2lm and 2mm- 

8.4 Transitional channel flow 

A known drawback of the traditional eddy viscosity closure for LES of transitioned 
flows is that it is overly dissipative, possibly eliminating instabilities altogether 
(Piomelli Zang, 1990). The dynamic model, on the other hand, yields essentially 
zero eddy viscosity if the resolved part of the flow is not turbulent. Instabilities 
are thus allowed to grow initially in a realistic fashion, as shown in simulations of 
transitional channel flow using the dynamic model, with planar averaging (Germano 
et a/., 1991), Once the non-linear breakdown phase is reached, the SGS model must 
become active in order to prevent excessive growth of turbulent kinetic energy, wall 
shear-stresses, etc. In the Lagrangian model, the variable 2cm must be initialized 
to zero everywhere in the laminar region. As turbulence is generated, this variable 
(and therefore the eddy viscosity) will rise from zero. The rate at which 2cm rises 
from zero is controlled in part by the memory time scale. If the memory time 
scale, T, is too long, the rise in eddy viscosity may occur too late in the transition 
process. In order to investigate this potential problem, we have performed an LES of 
transitional channel flow. In this section we attempt to ascertain if the Lagrangian 
model as proposed here (with the time scale given by Eq. 15) is able to (i) allow for 
initial instabilities to grow in a realistic fashion, and then to (ii) sufficiently damp 
the turbulence at the appropriate time. 

The transition channel case is identical to that of Zang et al. (1990), Piomelli & 
Zang (1991) and Germano et al. (1991). The initial (laminar) centerline Reynolds 
number is 8,000. The initial condition consists of a parabolic mean flow plus a 2-D 
Tollmien-Schlichting wave of 2% amplitude and a pair of 3-D Tollmien-Schlichting 
waves of 0.02% amplitude. The streamwise wavenumber for both the 2-D and 3-D 
modes is 1.0, whereas the spanwise wavenumber for the 3-D modes is ±1.5. See 
Zang et ai (1990) for more details on the initial conditions. The dimensions of the 
computational domain (streamwise, wall- normal, and spanwise) are 27r x 2 x 47r/3 
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FIGURE 16. Time history of wall-shear stress from the transitional channel sim- 
ulation. : Lagrangian model LES; : plane- averaged dynamic model 

LES; • : DNS of Zang et al. (1990). 



FIGURE 17. Streamwise velocity fluctuation profiles from the transitional channel 
flow simulation. Symbols; DNS of Zang et al. (1990). • : t = 176; * : t = 200; 

m : t = 220; : Lagrangian model LES; : plane- averaged dynamic 

model LES. 
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FIGURE 18. Reynolds shear stress profiles from the transitional channel flow 
simulation. Symbols: DNS of Zang et al (1990). • : t = 176; * : t = 200; 

m : t — 220; : Lagrangian model LES; : plane- averaged dynamic 

model LES. 

(in units of 6). The term Tim is initialized to 10 -14 (instead of to zero) in order 
to allow the first-order Euler scheme (explicit in T n ) to move Tim away from zero 
once the source term L : M becomes non-zero. 

The calculation is started on a 16 x 65 x 16 mesh. As the transition process 
proceeds, the solution is interpolated onto increasingly finer meshes. The timings 
of the remeshings are determined by monitoring the energy content in the highest 
resolved frequencies in the streamwise and spanwise directions. The remeshing 
procedure was found to introduce a complication in the Lagrangian SGS model. 
Refining the mesh while holding the test-to-grid filter ratio fixed results in different 
values of L : M and M : M. Because of its memory, the Lagrangian model requires 
a finite amount of time to adjust to the sudden changes in L : M and M : M (about 
A t = 5, or 100 timesteps). In order to minimize this recovery time, the remeshing is 
performed with values of Tim and Tmm rescaled so that their plane- averaged values 
are equal to those of the instantaneous L : M and M : M, respectively. Early in 
the transition process the SGS dissipation is minuscule and errors associated with 
remeshing probably have a negligible effect. However, the flow may be more sensitive 
to remeshing at later times when the SGS dissipation is not negligible. 

The 16 x 65 x 16 mesh is used until t = 145 (in units of initial centerline velocity 
U c and <5), when the grid is remeshed to 24 x 65 x 24. The run is then continued to 
t = 176 on both 24 x 65 x 24 and 32 x 65 x 32 meshes (with little notable difference). 
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The field is then remeshed to 32 x 65 x 48. Another remeshing to 48 x 65 x 64 is 
performed at t = 200, and the simulation is then run without further remeshing to 
t = 280. 

Figure 16 shows the time- history of the wall-shear stress compared with the DNS 
of Zang ei al. (1990). Results from the plane- averaged dynamic model are included 
in this figure. The Lagrangian model is in good agreement with the DNS results 
up to t = 210. Then, the wall-shear stress slightly overshoots the peak after which 
it settles to a plateau, near the DNS value. The plane-averaged dynamic model 
results axe similar, with the exception that the peak shear stress is underpredicted. 
Streamwise velocity fluctuations from the Lagrangian and plane-averaged models 
at times t = 176, 200, and 220 are compared with the (filtered) DNS data in Fig. 
17. Overall the agreement is quite good, and at t = 176 it is excellent. At this time 
the Lagrangian and plane-averaged results are indistinguishable. Reynolds shear 
stresses are shown in Fig. 18. Very good agreement is obtained at t = 176, whereas 
some differences exist at t = 200 and t = 220. 

Overall these results show that the Lagrangian model is capable of simulating 
transition. The eddy viscosity does rise from zero with a delay which is small 
enough so that turbulence is sufficiently damped after the rapid growth of kinetic 
energy during transition. 

4, Summary and conclusions 

A new version of the dynamic model has been tested in conjunction with the 
Smagorinsky closure. The model involves averaging the Germano identity for some 
time along fluid pathlines rather than over directions of statistical homogeneity, 
as was the practice in previous applications of the dynamic model. The present 
model is not restricted to flows with such special directions and should be readily 
applicable to complex-geometry, unsteady flows. We have shown that if an expo- 
nential memory is employed, the required averages can be obtained by solving a 
pair of relaxat ion- transport equations. In order to allow for the implementation of 
this model with minimal computational complications, we proposed to discretize 
the total derivatives that enter in these equations using a first-order expression in 
time, coupled with linear spatial interpolation to find the values required at the 
‘upstream’ locations. The resulting formulation (embodied in Eqs. (17) and (18)) 
is very simple to implement. 

Basic properties of the model were studied in DNS and LES of forced isotropic 
turbulence. The effect of the Lagrangian averaging on the pdfs of various quantities 
involved in the modeling were identified. It was also shown that the model preserves 
enough spatial locality to be influenced by vortical structures (‘fat worms’) that were 
identified in the LES. 

Applications of LES to isotropic turbulence and fully developed and transitional 
channel flow has shown that the model performs well and should be readily appli- 
cable to more complex flows. 

On a final note, we recognize that the Lagrangian dynamic model contains some 
arbitrary elements. In particular, an adjustable memory time scale T is involved. 
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This fact is unfortunate since it appears to conflict with the dynamic model philos- 
ophy of dispensing with adjustable parameters in favor of determining the subgrid- 
scale stress solely from information contained in the resolved velocity field. However, 
for any implementation which uses averaging, there is a similar ambiguity in choos- 
ing the domain over which Germano’s identity is to be enforced. Schemes that make 
use of spatial averaging often average over all homogeneous directions although a 
smaller subspace may be sufficient to insure the stability of the model. Choosing a 
particular value of the averaging time scale in the Lagrangian model is analogous 
to choosing a particular region in space over which to average. When viewed in 
this way, the Lagrangian model actually has an advantage over spatially-averaged 
variants in that it is designed to average over the minimum time necessary to insure 
stability. This feature allows the model to retain the maximum amount of spatial 
and temporal variability while remaining stable. 
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